{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local Homology\n",
    "\n",
    "In this tutorial we use tools from Topological Data Analysis (TDA) to study the local geometry of point clouds.\n",
    "\n",
    "If you are looking at a static version of this notebook and would like to run its contents, head over to [GitHub](https://github.com/giotto-ai/giotto-tda/blob/master/examples/local_homology.ipynb) and download the source.\n",
    "\n",
    "## See also\n",
    "\n",
    "- [Topological feature extraction using VietorisRipsPersistence and PersistenceEntropy](https://giotto-ai.github.io/gtda-docs/latest/notebooks/vietoris_rips_quickstart.html) for a quick introduction to general topological feature extraction in ``giotto-tda``.\n",
    "- [Local Homology NLP Use Cases: Unsupervised text disambiguation](https://giotto-ai.github.io/gtda-docs/latest/notebooks/local_hom_NLP_disambiguation.html), in which the ideas presented here are applied to natural language data!\n",
    "\n",
    "**License: AGPLv3**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from gtda.plotting import plot_point_cloud\n",
    "from gtda.local_homology import KNeighborsLocalVietorisRips, RadiusLocalVietorisRips"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dimensionality\n",
    "\n",
    "## What is \"dimensionality\", and how can we formalize it? \n",
    "\n",
    "In this tutorial we will learn how to algorithmically distinguish between possible *local geometries* around data points. Consider, for instance, the local geometry around the red dot in the three shapes below:\n",
    "\n",
    "![locally_different](images/locally_different.png)\n",
    "\n",
    "## How?\n",
    "\n",
    "We propose an answer to this question based on the notion of **local homology**. We estimate the \"local dimension\" around each point by studying the topological properties of the data points' neighbourhoods. We start by studying three shapes which differ by their dimensions, namely the line, the plane, and the cube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "n = 10\n",
    "line_1d = np.linspace(0, 1, num=n, endpoint=False)\n",
    "\n",
    "line = np.c_[np.zeros((n, 2)), line_1d]\n",
    "plane = np.array([[0, y, z] for y in line_1d for z in line_1d])\n",
    "cube = np.array([[x, y, z] for x in line_1d for y in line_1d for z in line_1d])\n",
    "\n",
    "line_plane_cube = np.vstack((line, plane + [0.5, -0.5, 0], cube + [1, -0.5, 0]))\n",
    "plot_point_cloud(line_plane_cube)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## The method\n",
    "\n",
    "As explained in [Topological feature extraction using VietorisRipsPersistence and PersistenceEntropy](https://giotto-ai.github.io/gtda-docs/latest/notebooks/vietoris_rips_quickstart.html), \"persistence diagrams\" are a common and useful way to store information about the topology and geometry of data. Their content can be summarized and made interpretable even further by means of a variety of featurization methods.\n",
    "\n",
    "In this notebook, we show how to compute and featurize persistence diagrams for *local neighbourhoods* of all vectors in a dataset (e.g. in a point cloud), in such a way that each point's features capture the dataset's dimensionality around that point. We do this in steps:\n",
    "1. create \"local point clouds\" around our data points and compute their persistence diagrams in a few chosen dimensions (more on this in \"How does it work?\" below);\n",
    "2. featurize the \"local persistence diagrams\" using a modified version of the *persistence entropy* computed by the ``PersistenceEntropy`` transformer;\n",
    "3. plot the original point cloud with its points coloured according to the value of the dimension-$i$ feature.\n",
    "\n",
    "Let us see this sequence of steps in action in the case of the point cloud representing a line. The colours will correspond to the value of the dimension-1 local topological feature. The local persistence diagrams are computed here using the ``KNeighborsLocalVietorisRips`` transformer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gtda.diagrams import PersistenceEntropy\n",
    "from sklearn.preprocessing import FunctionTransformer\n",
    "from gtda.pipeline import make_pipeline\n",
    "\n",
    "## 1 - Initialize a local homology transformer. \n",
    "# The notion of neighborhood here is determined by the number of closest neighbors to consider\n",
    "n_neighbors = (1, 7)\n",
    "homology_dimensions = (1, 2)\n",
    "kn_lh = KNeighborsLocalVietorisRips(n_neighbors=n_neighbors,\n",
    "                                    homology_dimensions=homology_dimensions)\n",
    "\n",
    "## 2 - Define a featurization method for persistence diagrams.\n",
    "mod_pe =  make_pipeline(PersistenceEntropy(),\n",
    "                        FunctionTransformer(func=lambda X: 2 ** X))\n",
    "\n",
    "# Summarize the steps above in a single end-to-end scikit-learn pipeline.\n",
    "# Then, apply the pipeline\n",
    "pipe = make_pipeline(kn_lh,\n",
    "                     mod_pe)\n",
    "loc_dim_features = pipe.fit_transform(line)\n",
    "\n",
    "## 3 - Plot the point cloud, colouring according to the value of the dimension-1 feature\n",
    "# A simple helper function\n",
    "def plot_coloured_cloud(point_cloud, colors):\n",
    "    \"\"\"This is a custom plotting function, that plots the point cloud and colours \n",
    "     the points according to a vector of colours.\"\"\"\n",
    "    plotly_params = {\"trace\" : {\"marker\": {\"color\": colors,\n",
    "                                           \"size\": 6,\n",
    "                                           \"showscale\": True},\n",
    "                                \"text\": [f\"color: {v}\" for v in colors]}}\n",
    "    return plot_point_cloud(point_cloud, plotly_params=plotly_params)\n",
    "\n",
    "# Find the column index for homology dimension 1 in `loc_dim_features`\n",
    "dimension = 1\n",
    "dim_index = homology_dimensions.index(dimension)\n",
    "colors = loc_dim_features[:, dim_index]\n",
    "\n",
    "# Plot the point cloud using the results in local homology dimension 1\n",
    "plot_coloured_cloud(line, colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each non-boundary point, the algorithm only detects one local feature, in dimension 1!\n",
    "\n",
    "Now let's try with the plane point cloud, colouring according to the value of the dimension-2 local topological feature.  The local persistence diagrams are computed here using the ``RadiusLocalVietorisRips`` transformer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plane coloured with two-dimensional features.\n",
    "# The notion of neighborhood here is determined by radii.\n",
    "radii = (0.001, 0.3)\n",
    "homology_dimensions = (1, 2, 3)\n",
    "r_lh = RadiusLocalVietorisRips(radii=radii,\n",
    "                               homology_dimensions=homology_dimensions)\n",
    "\n",
    "# Construct and apply the pipeline\n",
    "pipe = make_pipeline(r_lh,\n",
    "                     mod_pe)\n",
    "loc_dim_features = pipe.fit_transform(plane)\n",
    "\n",
    "# Find the column index for homology dimension 1 in `loc_dim_features`\n",
    "dimension = 2\n",
    "dim_index = homology_dimensions.index(dimension)\n",
    "colors = loc_dim_features[:, dim_index]\n",
    "\n",
    "# Plot the point cloud using the results in local homology dimension 1\n",
    "plot_coloured_cloud(plane, colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As for the line with its one-dimensional local features, the *two-dimensional* local features of the plane are equal to 1 everywhere except at the borders/boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-homogeneous example\n",
    "\n",
    "In the following example, we  consider a point cloud with varying dimensionality. We plot the three shapes next to each other. We encourage you to play around with the dimension parameter, to compare the local structure of points in the line, plane and cube!\n",
    "\n",
    "One sees that one-dimensional features are observed in all three shapes, but two-dimensional features only exist in the plane and three dimensional features only exist in the cube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we set the `collapse_edges` parameter to True, which often improves performance in higher homology dimensions.\n",
    "# We also set the `n_jobs` parameter to -1, to allow the algorithm to distribute across CPU cores!\n",
    "radii = (0.01, 0.25)\n",
    "homology_dimensions = (1, 2, 3)\n",
    "collapse_edges = True\n",
    "n_jobs = -1\n",
    "r_lh = RadiusLocalVietorisRips(radii=radii,\n",
    "                               homology_dimensions=homology_dimensions,\n",
    "                               collapse_edges=collapse_edges,\n",
    "                               n_jobs=n_jobs)\n",
    "\n",
    "# Construct and apply the pipeline\n",
    "pipe = make_pipeline(r_lh,\n",
    "                     mod_pe)\n",
    "loc_dim_features = pipe.fit_transform(line_plane_cube)\n",
    "\n",
    "# Find the column index for homology dimension 1 in `loc_dim_features`.\n",
    "# We encourage you to change the dimension parameter\n",
    "dimension = 1\n",
    "dim_index = homology_dimensions.index(dimension)\n",
    "colors = loc_dim_features[:, dim_index]\n",
    "\n",
    "# Plot the point cloud using the results in the chosen local homology dimension\n",
    "plot_coloured_cloud(line_plane_cube, colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \"Non-manifold\" points|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our last example is inherently one-dimensional, but the point cloud still has some interesting local structures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First generate the point cloud\n",
    "line1 = [[x, 0, 0] for x in np.arange(-1, 1, 1/10)]\n",
    "line2 = [[0, y, 0] for y in np.arange(-1, 1, 1/10)]\n",
    "line3 = [[x, 1, 0] for x in np.arange(-1, 1, 1/10)]\n",
    "lines = np.array(line1 + line2 + line3)\n",
    "\n",
    "# Plot the point cloud\n",
    "plot_point_cloud(lines)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, there are 4 different types of points:\n",
    "- 'edge points', with no local structure;\n",
    "- 'line points', which locally look like one-dimensional lines;\n",
    "- 'intersection points' which locally look like an X;\n",
    "- 'intersection points' which locally look like a T.\n",
    "\n",
    "The next block looks at the local homology features!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the radius local homology transformer.\n",
    "radii = (0.21, 0.31)\n",
    "homology_dimensions = (0, 1)\n",
    "n_jobs = -1\n",
    "r_lh = RadiusLocalVietorisRips(radii=radii,\n",
    "                               homology_dimensions=homology_dimensions,\n",
    "                               n_jobs=n_jobs)\n",
    "\n",
    "# Construct and apply the pipeline\n",
    "pipe = make_pipeline(r_lh,\n",
    "                     mod_pe)\n",
    "loc_dim_features = pipe.fit_transform(lines)\n",
    "\n",
    "# Find the column index for homology dimension 1 in `loc_dim_features`.\n",
    "dimension = 1\n",
    "dim_index = homology_dimensions.index(dimension)\n",
    "colors = loc_dim_features[:, dim_index]\n",
    "\n",
    "# Plot the point cloud using the results in the chosen local homology dimension\n",
    "plot_coloured_cloud(lines, colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the one-dimensional feature counts the branches at a given point!\n",
    "\n",
    "We encourage you to change the values of the ``radii`` parameter. If you are curious of what these parameters are, please scroll down to the 'How does it work?' section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How does it work?\n",
    "\n",
    "As the mathematics behind this algorithm is quite involved, we give an explanation for one-dimensional examples. We also give a picture for the process in a 2-dimensional example.\n",
    "\n",
    "A central theme in topology is the study of properties that are invariant under deformation. A popular one-dimensional invariant is the number of 'loops' in the space. For a graph, the number of 'loops' is the number of cycles of that graph. In the case of connected graphs this is equal to $1 - \\chi(G)$, where $\\chi(G) = |V| - |E|$ is the *Euler characteristic* of the graph, and $|V|$ and $|E|$ are the number of vertices and edges respectively. \n",
    "\n",
    "Let us consider *trees*, which are a particular kind of graph. Although all trees have Euler characteristic equal to $1$, they can have very different local structure. By computing this local structure, we can discriminate between topologically distinct trees (or graphs, or higher-dimensional objects). In order to compute the local structure around each data point with standard algorithms for point clouds or dissimilarity matrices, such as the highly-optimised routines used by ``VietorisRipsPersistence``, we use a few tricks!\n",
    "\n",
    "The first step is to discard points that are far away from the point under consideration, since they do not affect the local geometry around it. We thus only keep points that are \"not far\" from our point, according to a certain criterion. For the ``RadiusLocalVietorisRips`` transformer, the criterion is that the distance from the chosen point should be less than or equal to the second entry in the ``radii`` pair. For the ``KNeighborsVietorisRips`` transformer, we only include the $k$-nearest neighbors to our chosen point, where $k$ is the second entry in the ``n_neighbors`` pair.\n",
    "\n",
    "The second step is the one making loops reflect dimensionality. From a line, how can we obtain a loop? We stick the two ends of the line together! This trick can be thought of as \"coning off\" the endpoints: adding an extra point and linking some points (the furthest points) to this extra point. In the case of the ``RadiusLocalVietorisRips`` transformer, our algorithm cones off the points that lie in between the first and second entries in ``radii``, as depicted below. For the ``KNeighborsVietorisRips`` transformer, we cone off in an analogous way according to the ``n_neighbors`` pair.\n",
    "\n",
    "So what is the Euler characteristic of the result of applying this procedure to an \"X shape\" graph, which could be the local neighborhood around a point in a larger graph? The result after \"coning off\" has two vertices (at the intersection points) and 4 edges, so the Euler characteristic is $-2$! This reflects the complexity at the point that we are considering. (Notice also that a loop has Euler characteristic $0$, and a triple intersection will give Euler characteristic $-1$.)\n",
    "\n",
    "As a 2-dimensional example of this procedure, consider a disk (again, you can think of this as the local neighbourhood around a chosen point). By coning off the outer part of this disk, we end up with a sphere. And a sphere's only non-trivial topological features are 2-dimensional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![coning_off](images/coning_off.png)\n",
    "<div style=\"text-align: left\">\n",
    "   <p style=\"text-align: center;\">  Coning off a cross and a disk.</p>\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
